#!/usr/bin/python
# GM_LumFunc.py
'''Create realistic populations of stars, up to a constant.'''

import numpy

VariableNames = ['R', 'R0', 'hR', 'r', 'scale', 'galactocentric', 'ScaleList', 'LumFuncList', 'delMvList', 'z', 'hz']

def lumfunc(simpars, phypars):
    for i in simpars:
        cmd = "%s = simpars['%s']" % (i,i)
        exec cmd
    for i in phypars:
        cmd = "%s = phypars['%s']" % (i,i)
        exec cmd
    # R dependence of density, including a cutoff at r>scale to prevent assymetry about galactic center
    Rcorr = numpy.exp( -(R - R0)/hR )
    if galactocentric == True:
        rmask = numpy.greater(scale, r)
        Rcorr = Rcorr*rmask
    elif galactocentric == False:
        pass
    else:
        print "You're a cucumber!  ...No seriously, check the value you entered for 'galactocentric'.  It should be boolean."
    distriblist = []
    del rmask
    # z dependence of density
    for i in range(0, len(MvList)):
        hz = ScaleList[i]
        rho0Mv = LumFuncList[i]
        binwidth = delMvList[i]
        zcorr = numpy.exp(-z/hz)
        thisdistrib = binwidth*rho0Mv*numpy.multiply(Rcorr, zcorr[:,numpy.newaxis,:])
        distriblist.append(thisdistrib)
    del Rcorr, zcorr
    return distriblist
